### Sweden WEP R and R 

library(tidyverse)

library(lubridate)
library(broom)
library(gganimate)
library(transformr)
library(gtools)
library(jtools)
library(huxtable)
library(estimatr)
library(texreg)
library(modelsummary)


# Using up to August

se_local <- read.csv("Sweden_mobility.csv")



Sweden_2018 <- read_csv("Sweden_2018.csv")

Sweden_2018 <- Sweden_2018 %>% 
  mutate(municipality = str_sub(municipality, 6L, -1L))

Sweden_2018<-Sweden_2018 %>% 
  mutate(sub_region_2 = paste(municipality, "Municipality"))

se_local<-se_local %>% 
  left_join(Sweden_2018, by  = "sub_region_2")

se_local <- se_local %>% 
  filter(sub_region_2!="")

se_local  <- se_local %>% 
  mutate(date = as.Date(date)) %>% 
  rename (grocery = grocery_and_pharmacy_percent_change_from_baseline,
          parks = parks_percent_change_from_baseline,
          residential = residential_percent_change_from_baseline,
          retail_rec = retail_and_recreation_percent_change_from_baseline,
          transit = transit_stations_percent_change_from_baseline,
          workplace = workplaces_percent_change_from_baseline,
          location = sub_region_1) %>% 
  mutate(id_no = as.numeric(id))

se_covid <- read_csv("veckodata_kommun.csv")

se_covid <- se_covid %>% 
  filter(veckonummer<32) %>% 
  filter(KnKod!="Okänd") %>% 
  mutate(antal_fall_per10000_inv = case_when(is.na(antal_fall_per10000_inv) ~ 0, 
                                                TRUE ~ antal_fall_per10000_inv),
         date = ((veckonummer*7)+as.numeric(as.Date("2019-12-29"))),
         date = as.Date(date, origin = "1970-01-01"),
         id = KnKod,
         id_no = as.numeric(KnKod))

date_data <- se_local %>% 
  select(c(id, date))

se_covid <-  se_covid %>% 
  full_join(date_data, by = c("id", "date"))


se_covid_merge  <- se_covid %>% 
  group_by(id, date) %>% 
  summarise(week_cases = mean(antal_fall_per10000_inv, na.rm=T)) %>% 
  mutate(id = str_sub(id, -4, -1))

library(zoo)

se_local <-  se_local %>% 
  left_join(se_covid_merge, by=c("id", "date"))%>% 
  group_by(sub_region_2) %>% 
  arrange(sub_region_2, date) %>% 
  mutate(week_cases_i  = na.approx(week_cases, maxgap = 6, rule = 2))



### 

# se_full_year <- read.csv( "https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv")
# 
# 
# se_full_year <- se_full_year %>% 
#   filter(country_region_code=="SE")
# 
# se_full_year<-se_full_year %>% 
#   left_join(Sweden_2018, by  = "sub_region_2")
# 
# se_full_year <- se_full_year %>% 
#   filter(sub_region_2!="")
# 
# se_full_year  <- se_full_year %>% 
#   mutate(date = as.Date(date)) %>% 
#   rename (grocery = grocery_and_pharmacy_percent_change_from_baseline,
#           parks = parks_percent_change_from_baseline,
#           residential = residential_percent_change_from_baseline,
#           retail_rec = retail_and_recreation_percent_change_from_baseline,
#           transit = transit_stations_percent_change_from_baseline,
#           workplace = workplaces_percent_change_from_baseline,
#           location = sub_region_1)

# Sweden Graphs

weekend_df_se <- data.frame(date = unique(se_local$date)) %>%
  mutate(day = weekdays(date)) %>%
  mutate(weekend= case_when(
    day %in% c("Saturday", "Sunday") ~ 1,
    date==as.Date("2020-04-10") ~ 1,
    date==as.Date("2020-04-13") ~ 1,
    date==as.Date("2020-05-01") ~ 1,
    date==as.Date("2020-05-21") ~ 1,
    date==as.Date("2020-05-22") ~ 1,
    date==as.Date("2020-06-19") ~ 1,
    TRUE ~ 0))

coefs_map_workplace_se <- map_df(unique(se_local$date), function(date_i) {
  filtered_data <- filter(se_local, date == date_i) %>%
    lm(data =., workplace ~ sdshare+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_workplace_se <- coefs_map_workplace_se %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("sdshare", "medianincome", "employmentrate", "populationdensity", "postsecond", "price_villas")) %>% 
  mutate(term = case_when(term=="sdshare" ~ "Sweden Democrats Vote",
                          term=="medianincome" ~ "Median Income",
                          term=="employmentrate" ~ "Employment",
                          term=="populationdensity" ~ "Density",
                          term == "postsecond" ~ "Post Secondary Education",
                          term == "price_villas" ~ "House Prices"))


coefs_map_workplace_se_cases <- map_df(unique(se_local$date), function(date_i) {
  filtered_data <- filter(se_local, date == date_i) %>%
    lm(data =., workplace ~ sdshare+week_cases_i+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_workplace_se_cases <- coefs_map_workplace_se_cases %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("sdshare", "medianincome", "week_cases_i", "populationdensity", "postsecond", "price_villas")) %>% 
  mutate(term = case_when(term=="sdshare" ~ "Sweden Democrats Vote",
                          term=="medianincome" ~ "Median Income",
                          term=="week_cases_i" ~ "Weekly Cases",
                          term=="populationdensity" ~ "Density",
                          term == "postsecond" ~ "Post Secondary Education",
                          term == "price_villas" ~ "House Prices"))


se_local %>% 
  ggplot(aes(x = date, y = workplace, group=sub_region_2, alpha = 0.2))+
  geom_jitter(aes(color=as.factor(sdshare>20.5)))+
  geom_vline(xintercept =as.Date("2020-03-25"), linetype="dashed", color = "red")+
  scale_color_manual(values=c('yellow','blue'))+
  ylab("Change in Workplace Activity from Baseline")+xlab("Date")+
  theme_classic()+
  theme(legend.position = "none")

ggsave("se_dragon.pdf", width = 8, height = 5)
ggsave("se_dragon.jpg", width = 8, height = 5)

coefs_map_workplace_se %>% 
  left_join(weekend_df_se) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("se_coefs.pdf", width = 8, height = 5)
ggsave("se_coefs.png", width = 8, height = 5)

coefs_map_workplace_se_cases %>% 
  left_join(weekend_df_se) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("se_coefs_cases.pdf", width = 8, height = 5)
ggsave("se_coefs_cases.png", width = 8, height = 5)




# Sweden regressions ------------------------------------------------------


se_local <- se_local %>% 
  mutate(after_march_16 = as.numeric(date>as.Date("2020-03-16")),
         after_march_25 = as.numeric(date>as.Date("2020-03-25")),
         #  after_sep_13 = as.numeric(date>as.Date("2020-09-13")),
         #   after_nov_4 = as.numeric(date>as.Date("2020-11-04")),
         #  after_dec_1 = as.numeric(date>as.Date("2020-12-01"))
  ) %>% 
  group_by(location) %>% 
  mutate(lag_workplace = lag(workplace),
         l2_workplace = lag(workplace, 2),
         l3_workplace = lag(workplace, 3),
         l4_workplace = lag(workplace, 4),
         l5_workplace = lag(workplace, 5),
         l6_workplace = lag(workplace, 6)) %>% 
  mutate(day = weekdays(date)) %>%
  mutate(weekend= case_when(
    day %in% c("Saturday", "Sunday") ~ 1,
    date==as.Date("2020-04-10") ~ 1,
    date==as.Date("2020-04-13") ~ 1,
    date==as.Date("2020-05-01") ~ 1,
    date==as.Date("2020-05-21") ~ 1,
    date==as.Date("2020-05-22") ~ 1,
    date==as.Date("2020-06-19") ~ 1,
    TRUE ~ 0)) %>% 
  ungroup()


se_weekdays <- se_local %>% 
  filter(weekend==0)



s3<-lm_robust(data = se_weekdays, workplace ~ lag_workplace+  sdshare*after_march_25+
              medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas)

s4<-lm_robust(data = se_weekdays, workplace ~ lag_workplace+  sdshare*after_march_25+
                medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
              fixed_effects = location)

s5<-lm_robust(data = se_weekdays, workplace ~ lag_workplace+  sdshare*after_march_25+
                medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
              fixed_effects = location+date)

r1 <- lm_robust(data = se_weekdays, workplace ~ lag_workplace+
                  (sdshare)*after_march_25,
                #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                fixed_effects = location+date)

r1a <- lm_robust(data = se_weekdays, workplace ~ lag_workplace+
                  (sdshare+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas)*after_march_25,
                #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                fixed_effects = location+date)

r2 <- lm_robust(data = se_weekdays, workplace ~ lag_workplace+
                  (sdshare)*after_march_25,
                #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                fixed_effects = sub_region_2+date)

r2a <- lm_robust(data = se_weekdays, workplace ~ lag_workplace+
                  (sdshare+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas)*after_march_25,
                #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                fixed_effects = sub_region_2+date)


coefs_s3 <- list( "lag_workplace" = "Lag Workplace", 
                  "sdshare" = "Sweden Democrats Vote" , "medianincome" = "Median Income",
                  "employmentrate" = "Employment Rate" , "populationdensity" = "Population Density",
                  "primary" = "% Primary Educated", "uppersecondary" = "% Upper Secondary Educated" ,
                  "price_villas" = "House Prices" ,  "after_march_25" = "After March 25th" ,
                  "sdshare:after_march_25" = "SD Share * After Mar 25")


# This table has Model 1 
texreg(list(s4, s5,  r1a, r2a), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "County Effects"= c("Y", "Y",  "Y",  "N"), "Date Effects" = c( "N", "Y",  "Y",  "Y"),
                               "Muncipality Effects" = c("N", "N",   "N",  "Y"), "All Interactions" = c("N", "N",  "Y",  "Y"))
       , file="Sweden_RandR.tex", caption = "Workplace Activity in Sweden (no weekends)", 
       caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_se", booktabs = TRUE)


se_weekdays_l <- se_weekdays %>% filter(medianincome<=25.4)

se_weekdays_h <- se_weekdays %>% filter(medianincome>25.4)


r2a_l <- lm_robust(data = se_weekdays_l, workplace ~ lag_workplace+
                   (sdshare+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas)*after_march_25,
                 #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                 fixed_effects = sub_region_2+date)

r2a_h <- lm_robust(data = se_weekdays_h, workplace ~ lag_workplace+
                     (sdshare+medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas)*after_march_25,
                   #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                   fixed_effects = sub_region_2+date)

s5l<-lm_robust(data = se_weekdays_l, workplace ~ lag_workplace+  sdshare*after_march_25+
                 medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
               fixed_effects = location+date)

s5h<-lm_robust(data = se_weekdays_h, workplace ~ lag_workplace+  sdshare*after_march_25+
                 medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
               fixed_effects = location+date)


summary(r2a_l)

summary(r2a_h)

coefs_s3 <- list( "lag_workplace" = "Lag Workplace", 
                  "sdshare" = "Sweden Democrats Vote" , "medianincome" = "Median Income",
                  "employmentrate" = "Employment Rate" , "populationdensity" = "Population Density",
                  "primary" = "% Primary Educated", "uppersecondary" = "% Upper Secondary Educated" ,
                  "price_villas" = "House Prices" ,  "after_march_25" = "After March 25th" ,
                  "sdshare:after_march_25" = "SD Share * After Mar 25")


# This table has Model 1 
texreg(list(s5l, s5h, r2a_l, r2a_h), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "Income Level" = c("Low", "High", "Low", "High"), "Date Effects" = c( "Y", "Y", "Y",  "Y"),
                               "Muncipality Effects" = c( "N", "N","Y",  "Y"), "All Interactions" = c( "N", "N", "Y",  "Y"))
       , file="Sweden_RandR_split.tex", caption = "Workplace Activity in Sweden (no weekends): Split by Income", 
       caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_se_split", booktabs = F)


s4l<-lm_robust(data = se_weekdays_l, workplace ~ lag_workplace+  sdshare*after_march_25+
                medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
              fixed_effects = location)

s4h<-lm_robust(data = se_weekdays_h, workplace ~ lag_workplace+  sdshare*after_march_25+
                 medianincome+employmentrate+populationdensity+primary+uppersecond+ postsecond+price_villas,
               fixed_effects = location)

summary(s4l)
summary(s4h)
